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Abstract 



We reconsider the theory of the linear response of non-equilibrium 
steady states to perturbations. We first show that by using a general 
functional decomposition for space-time dependent forcings, we can 
define elementary susceptibilities that allow to construct the response 
of the system to general perturbations. Starting from the definition of 
SRB measure, we then study the consequence of taking different sam- 
pling schemes for analysing the response of the system. We show that 
only a specific choice of the time horizon for evaluating the response of 
£N1 the system to a general time-dependent perturbation allows to obtain 

the formula first presented by Ruelle. We also discuss the special case 
of periodic perturbations, showing that when they are taken into con- 
l/") sideration the sampling can be fine-tuned to make the definition of the 

correct time horizon immaterial. Finally, we discuss the implications 
of our results in terms of strategies for analyzing the outputs of numer- 
ical experiments by providing a critical review of a formula proposed 
by Reick. 

c3 1 Introduction 



The study of how the properties of general non-equilibrium statistical 
mechanical systems change when considering a generic perturbation, 
usually related to variations either in the value of some internal param- 
eters or in the external forcing, is of great relevance, both in purely 
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mathematical terms and with regards to applications to the natural 
and social sciences. Whereas in quasi-equilibrium statistical mechan- 
ics it is possible to link the response of a system to perturbations to 
its unforced fluctuations thanks to the fluctuation-dissipation theorem 
[91 [28], in the general non-equilibrium case it is not possible to frame 
rigorously an equivalence between internal fluctuations and forcings. 
At a fundamental level, this is closely related to the fact that forced 
and dissipative systems feature a singular invariant measure. Whereas 
natural fluctuations of the system are restricted to the unstable mani- 
fold, because, by definition, asymptotically there is no dynamics along 
the stable manifold, perturbations will induce motions - of exponen- 
tially decaying amplitude - out of the attractor with probability one, 
as discussed in, e.g., [231 Ell US] • It is worth noting that Lorenz antic- 
ipated some of these ideas when studying the difference between free 
and forced variability of the climate system [T2] . This crucial difficulty 
inherent to out-of-equilibrium systems is lifted if the external perturba- 
tion is, rather artificially, everywhere tangent to the unstable manifold, 
or if the system includes some stochastic forcing, which smooths out 
the resulting the invariant measure [10]. 

Recently, Ruelle E01 EH E3] paved the way to the study of the re- 
sponse of general non-equilibrium systems to perturbations by present- 
ing rigorous results leading to the formulation of a response theory for 
Axiom A dynamical systems [TS], which possess a Sinai- Ruclle-Bowen 
invariant measure |26j . Given a measurable observable of the system, 
the change in its expectation value due to an e-perturbation in the 
flow (or in the map, in the case of discrete dynamics) can be written 
as a perturbative series of terms proportional to e", where each term of 
the series can be written as the expectation value of some well-defined 
observable over the unperturbed state. Ruelle's formula is identical to 
Kubo's classical formula [8] when a Hamiltonian system is considered 

Whereas Axiom A systems are mathematically non-generic, the 
applicability of the Ruelle theory to a variety of actual models is sup- 
ported by the so-called chaotic hypothesis [5] , which states that systems 
with many degrees of freedom behave as if they were Axiom A systems 
when macroscopic statistical properties are considered. The chaotic 
hypothesis has been interpreted as the natural extension of the classic 
ergodic hypothesis to non- Hamiltonian systems [5]. 

In the last decade great efforts have been directed at extending 
and clarifying the degree of applicability of the response theory for 
non-equilibrium systems along five main lines: 

• extension of the theory for more general classes of dynamical 
systems [HE]; 

• introduction of effective algorithms for computing the response 
in dynamical systems with many degrees of freedom [I], in order 
to support the numerical analyses pioneered by |18l |3] ; 

• investigation of the frequency-dependent response - the suscep- 
tibility - for the linear and nonlinear cases, with the ensuing in- 
troduction of a new theory of Kramers-Kronig relations and sum 
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rules for non-equilibrium systems [131 125] supported by numerical 
experiments [13] : 

• study of the response to external perturbations of non-equilib- 
rium systems undergoing stochastic dynamics [17] 127] ; 

• use of Ruelle's response theory to study the impact of adding 
stochastic forcing to otherwise deterministic systems [15] . 

In particular, the response theory seems especially promising for tack- 
ling notoriously complex problems such as those related to studying 
the response of geophysical systems to perturbations, which include 
the investigation of climate change; see discussions in [U [14] [16] . In 
particular, in [16] . it is discussed that by deriving from the linear sus- 
ceptibility the time-dependent Green function, it is possible to devise 
a strategy to compute climate change for a general observable and for 
a general time-dependent pattern of forcing. Recently, response theory 
is becoming of great interest also in social sciences such as economics 

When developing a response theory, there are two possible ways 
to frame the temporal impact of the additional perturbation to the 
dynamics. Either one considers the impact at a given time i of a 
perturbation affecting the system since a very distant past, or one 
considers the impact in the distant future of perturbations starting at 
the present time. When deriving the response formula, Ruelle takes the 
first approach and delivers the correct formula [20] . Taking a different 
point of view and considering the specific case of periodic perturbations 
- which, anyway, tell us the whole story about the response by linearity 
-, Reick [18] derives a formula that is well suited for analyzing the 
output of numerical experiments |14| 116] . 

Given the great relevance and increasing popularity in applications 
of the response theory introduced by Ruelle, in this paper, we recon- 
sider the theory of the linear response of non-equilibrium steady states 
to perturbations and try to bridge the theoretical derivations and the 
strategies for designing numerical experiments and analyzing efficiently 
their outputs. 

In Section 2, we study the relevance of the choice of the time horizon 
for evaluating the impact of the perturbation and we demonstrate by 
direct calculation that the Ruelle approach is the correct one. We 
clarify some of the assumptions implicitly considered in his derivations. 
We then discuss the special case of periodic perturbations, showing that 
using them as basis for a response theory greatly simplifies the formulas 
and the conditions under which the formulas are derived. In Section 
3, we discuss the implications of our results in terms of strategies for 
improving the quality of numerical simulations and of the analysis of 
their output signals and reconsider Reick's formula [18]. In Section 4 
we present our conclusions and perspectives for future work. 
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2 Linear Response Theory, revised 



2.1 Separable perturbations 

We study the linear response of a discrete dynamical system to general 
time-dependent perturbations. All calculations arc formal, in the sense 
that we neglect all higher orders in the perturbation without deriving 
an estimate for these terms and we assume that all sums converge in 
all senses necessary. 

The unperturbed dynamical system is given by 

xt+i = f(x t ) , 

with t £ Z, xt £ M, M being a smooth manifold and / : M — > M 
a differentiable map. For simplicity we consider a time-independent 
unperturbed dynamics, although the following can be extended to a 
time-dependent case in a straightforward manner. Moreover, the anal- 
ysis of the case of a continuous time flow x = f(x) is perfectly analo- 
gous to what is presented in the following and the main corresponding 
results will be mentioned in Appendix [Aj 

The dynamical system is perturbed by a time-dependent forcing 
X(t, x) as follows: 

x t+1 = f t+1 (x t ) := f(x t ) + X(t + 1, f(x t )) . (1) 

The effect of the perturbation on individual trajectories is in general 
difficult to describe. More can however be said about the statistical 
properties of the system. One can look at the expectation values of 
observables under invariant states of the dynamical system: 

P (A) ■= J p(dx)A{x) , 

where p(dx) is an invariant measure of the unperturbed dynamics i.e. 

P (Aof) = p(A). 

for any observable A. In general, a dynamical system can possess many 
invariant measures. The physically relevant measure for dynamical 
systems is the SRB measure |26j . This measure is physical in the sense 
that for a set of initial conditions of full Lebesgue measure the time 
averages lim^oo | 53fe=i A{f k (x)) converge to the expectation value 
under p. In other words, for any measure l(dx) that is absolutely 
continuous w.r.t. Lebesgue, we have that 

P(A) =\im -/t I l(dx)A(f k (x)). (2) 

°° k=l J 

We want to determine the linear response of expectation values 
under the SRB measure to perturbations of the dynamical system as 
in Eq. [T] We denote by 5tP the difference in the expectation value 
between the perturbed and unperturbed system at time T. In [2"5] . 
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Ruelle presents a formula for the linear response due to perturbations 
that are separable in time and space: 

X{t,x) = <f>{t) X {x). 

The leading order term of the expansion of Stp(A) in X is given by 

8 T p{A)^Y, G ^)4>{T-3), (3) 

with 

G A (j) = 0(j) J p(dx) X (x)D(A o f)(x) , (4) 

where 9 is the Heaviside function. Since 8tp{A) is expressed as a con- 
volution product of G a and 0, the Fourier transform of the response 
Suip(A) = J2tgz e lTul STp(A) is given by a product of the Fourier trans- 
form 4>{uj) of the time factor (f>(t) and a susceptibility function ka(w): 

5 u p(A) « k A (w)$(u) . (5) 



E ei3 ' u [ p(dx) X (x)D(Aof)(x). (6) 

Due to the causality of the response function G A (j) (i-e. Ga{]) = 0, 
j < 0)), the susceptibility ^(w) is analytic in the upper complex plane 
and satisfies Kramers-Kronig relations [HI [13] . 

2.2 General perturbations 

If the perturbation is of a more general nature (i.e. not separable), we 
can deduce a linear response formula from Eq. [5] solely based on lin- 
earity in the following way. Let <j> r (t) be a Schauder basis [11] of time- 
dependent functions and ip s {x) a Schauder basis of space-dependent 
functions. One can take for example the Fourier basis in time and 
a wavelet basis in space, or whatever basis may be suitable for the 
system at hand. The product functions 4> r {i)Tp s {x) then form a basis 
of the time and space dependent functions as a tensor product |24j . 
More concretely, we may for an appropriate sense of convergence as- 
sume that any function X(t, x) can be decomposed in the product basis 
4> r (t)ip s (x) with coefficients a r>s : 

X(t,x)= 2J a r . s (t) r (t)ip s {x) . 

r,s>0 



where 

k*>) = 
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Since each of the factors in this sum is separable, we can use Eq. [5] and 
the linearity of the response to get that the response is given by 

a r , s <f>r(u>)K St A(u) , (7) 

r,s>0 

where k s ^ a is the susceptibility function of observable A, corresponding 
to the forcing pattern given by tp s {x). Since the vectors ip s {x) consti- 
tute a basis, the functions k Sta are elementary linear susceptibilities 
that allow to construct the response of the system to any pattern of 
forcing. 

By inserting the expression of k Sj a(w) from Eq. [6] into Eq. [7j it 
is possible to deduce the frequency-dependent response of the system. 
It is expressed as an ensemble average of a dot product of Fourier 
transforms, namely the transforms of the perturbation term and of the 
linear tangent of the observable, G A (u,x): 

S u p(A) « ]Te^ [ p(dx)X(u,x)D(Aof)(x) 

= f p(dx)X(u),x)G A (uj,x), (8) 

with 

£(u,x) = J £ i X{T,x)e i " T 

= ^ a r , s (i>r(^)ips(x). (9) 

r,s>0 

Instead from Eqs.[3]|4]in the time domain 

Stp(A)^ f p(dx)Y,X(T-j,x)D(Aof)(x). (10) 

For the case of a periodic perturbation X(t + t,x) = X(t,x), where 
t G N, we get as linear response 

r oo 



p(dx)J2 X{T - n - m,T,x)D{Ao f n+mT ){x) 

n—1 m=0 

= I p(dx) X(T - n, x)G A<n (x) , (11) 

^ n=l 

with 

oo 

G A , n {x)=Y,D{Aop+™){x). 
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In order to elucidate some crucial aspects of the Ruelle's response 
theory, we now propose a direct derivation of the linear response to 
the perturbation X(t,x) by considering the history of the perturbed 
and unperturbed trajectory of the system and verify under which con- 



ditions we find agreement with Eqs. |8fl0| Our goal is to derive the 
leading order term of the expansion of Stp(A) with respect to X from 
first principle, i.e. without resorting to the Schauder decomposition as 
above. Such a derivation should of course arrive at the same results as 
those in Eqs. [MO 



2.2.1 Response at a moving time horizon 

We describe the perturbed measure Pt{A) such that the system is 
initialized at time T in an initial condition according to the measure I. 
We move the time horizon at which we observe forward and average the 
time-evolved measurements. The system is prepared and then observed 
while it is evolving over a sufficiently long time. The measure fix is 
time-dependent as the dynamics / is also time-dependent. Formally 
we take px to be the ergodic mean of the expectation values of A, 
starting at time T: 

p T (A) = lim i W l{dx)A{f T (x)) . (12) 

t^oot^J 

Here l(dx) is an initial measure that is absolutely continuous with 
respect to Lebesgue and fif, represents k iterations of the perturbed 
dynamics from time T to T + k: 

f T {x) = f T+k a...of T+l {x). (13) 

The difference in expectation values StP is the given by 

6 T p(A) = p T {A) - p(A). (14) 

Following the computation presented in [23] for the separable case, 
we can expand the perturbed dynamics / around the unperturbed 
dynamics /. We then try to rewrite the response of the perturbed 
system in terms of the SRB measure of the unperturbed system by 
finding an expression for A{fj,{x)) in terms of A(f k (x)). 

We can approximate up to first order in X the two time step future 
evolution by expanding around the unperturbed dynamics f 2 {x): 

xt+2 = fr+2 ° /t+i(£t) 

« f(x T ) + X(T + 1, f(x T )).Df(f(x T )) + X(T + 2, f(x T )) . 

For k time steps we similarly get: 

xr+k = fr+k 0...0 f T +i(x T ) 

k 

* f k i%T) +J2 X ( T + 3, P(x T )).(Df k -3)(P(x T )) . 
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Thus, we can approximate A(fo+k ° ■ ■ ■ fT+i(%)) to first order in X 
as follows: 

A(f T+k o...of T+1 (x))^A(f k (x)) 

+ A'(f k (x)) (^X(T + j,P(x)).(Df k -i)(f(x)) 
= A{f\x)) 

k 

+ J2x(T + j, f(x))D(A o f (x)) . 

3=1 

(15) 

The linear response of A is obtained by substituting Eq. [15] into 
Eq. [14} through Eq. [2] and Eq. [12] 

S TP (A) « Um - I l(dx)J2 X ( T + jJ j (x)M A °f k - j )(f j (x)) 

k=l J j=l 
t—1 t—i „ 

= t 1 ™ 7 EE / Kdx)X(T + j,f(x))D(Aof)(f(x)). 
Using that for i > t the expression is zero, we have 

5 TP (A)^Y.J [^\f^ryKdx)x{T + ^0^. 

(16) 



i>0 \ 3=1 



Note that it is not possible to rewrite the sum in j as the ergodic time 
mean of I due to the time dependence of the perturbation X(T + j, x). 
Therefore, surprisingly, Eq. [16] does not in general agree with Eq. [TU] 
In particular, by taking the limit on the right hand side, we obtain 
that the T-dependence disappears. Say we shift T to T — T' in the 
limit appearing in the above equation: 



lim \ y2{f*yi(dx)X{T-T' +j,x) 

3=1 
1 t-i-T' 

= lim- V (f*y'+ T 'l(dx)X(T + f,x) 

£— s-oo t * — ' 

j'=l-T> 

!™U £(/*)*' ((f*fKdx))x(T + j',x). 



t->i 



Taking the reasonable assumption that the result in Eq. [16] does not 
depend on the initial measure l(dx) (this cannot be obtained from the 
uniqueness of the SRB measure), the obtained response of the system 
is time-independent even if the forcing is time-dependent. 
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Let us compare the result contained in Eq. [16] with Eq. [3] in the 
special case of a time-independent perturbation X(t,x) = x( x )- Now 
Eq. [16] and Eq. [3] agree since Eq. |T6] simplifies to: 

S TP (A)^Y1 I p(dx) X (x)D(Aof )(x), 



j>0 



because lim^oo V*2j=i(/*)' J "'(^ a ')) = P(d x ), by the definition of the 
SRB measure. The formula given by Ruelle |21j is recovered, as can 
be seen by substituting tf>(t) = 1 into Eq. [3] 

However, already in the case of a time-periodic perturbation 

X(t,x) = X(t + r,x) 

there is no agreement between Eq. [16] and Eq. |10| In this case the sum 
over j appearing in Eq. [16] can be written as a double sum, one over 
k periods, indexed by m, and one over the r phases in each period, 
indexed by n: 



lim i y2(f*yi(dx)X(T + j,x) 

3=1 

1 k t 

= lim - Y,(n mT (f) n l(dx)X(T + n,x) 

m— 1 ri— 1 

i r 

-Y,Pn{dx)X(t + n,x), (17) 



T n=l 

with 

k 

p n (dx)= lim Tj^U^VTKdx)- 

m— 1 

Under the assumption that p n = p for all n G {l,...,r}, i.e. sub- 
sampling does not impact the unperturbed invariant measure, the re- 
sponse gives a similar result as the Ruelle formula, but with an av- 
eraged perturbation. Substituting Eq. [TT] into Eq. |16[ we obtain a 
formula of the form of Eq. [lOj with the difference that instead of the 
true forcing X(t,x) the averaged forcing 



- ^ X(t + n, x) 



T 

n=l 



appears. The disagreement is apparent, e.g. when one considers a per- 
turbation of the form X(t,x) = sin(— t)x(%), which obviously results 
in a zero response. This effect has a clear intuitive interpretation. The 
response at a given time depends mostly on the immediate past, hence 
if one does not keep fixed the horizon, one risks to average out the 
variability. The previous formula reflects this intuition. 



One way to obtain agreement with Formula 11 is to choose a specific 



sampling procedure. We sample with the same periodicity r of the 



9 



forcing, thus altering the definition of the response. We define the 
measures for the perturbed and unperturbed system as 

1 N 

p' T>p (dx) := lim ^Y,(fT +P ri(dx) 

iv— >oo 1\ L — ' 

fc=0 

1 N 

p'(^):= lim -£(/ fcT+I T«- (18) 

fc=0 



With this definition we obtain using Eq. 15 
5p' T ^A) = p' TtP {A) p'(A) 

N / N-l N-m 



^ J™, v E ( E E °( mT + n +p)) I (r T+n+v Y Kdx) 

~ >0 ° fc=o Vm=-1 i=l / J 

X(T + mT + n+ Pl x)D(A o f kr - mT - n )(x) . 

Using the periodicity of X one can obtain 

^T, P (^)«EE / UToo^ E (f n+mT+P ri(dx)e(mr + p + n)\ 

n=l i>l J \ m=-l / 

X(T + n+p,x) (D{Ao r- n )(x)) 

T-l 

^ / p{dx)X(T + p - n, x)G A ^ n (x) = S T+p p(A) . 

n J 



n=0 ' 



Hence by choosing the initial phase p at which we start sampling, we 
can obtain the response at this phase. This means that we only need 
to start one long simulation of / and / and d o s ummations of the 



differences (A o / — A o f){x) according to Eq. 18 at all phases p in 
one period to obtain the entire response to the periodic forcing. By 
applying a forcing that contains several frequencies, such as a block 
wave, we can extract the susceptibility at all present frequencies in one 
run by taking the Fourier transform of the response. 

Note that if we sample the signal with a periodicity r\ which is 
prime with respect to the period r of the forcing, we will obtain no 
p-dependence (with p, in this case, ranging from to 77 — 1) in the 
response. For all values of p we will obtain as a result the response to 
the time-averaged forcing. Therefore, the case of sampling at all time 
steps discussed above is just the special case given by 77 = 1, where we 
are basically considering the case of the Nyquist frequency. Instead, 
if t and 77 are not prime with respect to each other, the sampling 
procedure will be able to ascertain the p-dependence of the response 
of the system at the periodicity given by the common harmonic terms. 

If the periodicity of the forcing is not known, the above discussion 
tells us that by doing a sampling at larger and larger periods 77 and 
checking for each of those the phase-dependence of the response, it is 
possible to deduce the fundamental period of the forcing. If the proce- 
dure does not converge, we are facing a quasi-periodic or continuous- 
spectrum forcing for which this approach fails. 
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Therefore, this situation is unsatisfactory. Why do we only get the 
correct result for periodic perturbations and fine-tuning the sampling 
or by taking constant perturbations? 

2.2.2 Response at a fixed time horizon 

This paradox can be resolved by defining the time-dependent SRB 
measure in Eq. [12] using a different method of sampling. We now 
consider the time evolution in this definition to go from time T — k 
in the past up to the fixed time horizon T, so instead of Eq. |13| we 
have: 

f T = fa o . . . o / T _ fc . (19) 

Note that this approach does not use the reversed time dynamics but 
rather a different time perspective in which the final time is fixed as 
the current time and the perturbation starts in the remote past. 

The expansion to first order in X around the dynamics of / now 
becomes: 

XT = fa ° ■ ■ ■ ° fa-k+l{xT-k) 

fe-1 

« f\x T -k) + J2 X ( T - h f k - j (xT-k)).(Df)(f k - j (x T . k )) . 
3=0 

(20) 

Hence, the linear response of p(A) at time T is given by 

5 TP (A) *\im - t J2 I l(dx)J2x(T-j,f k -i(x))D(Aop)(f-j(x)) 

= t l™ \ £ / Y,(ryi(dx)x(T - h X )d{a p){ X ) . 



Note that in contrast to Eq. 16 the indices are such that the time aver- 
age of the measure and the perturbation are decoupled. This crucially 
depends on the choice of the sampling. This allows us to use the def- 
inition of the SRB measure in Eq. [2] and replace the time average in 
the limit by p: 

8 TP (A)^J2 I p(dx)X(T-j,x)D(Aof)( x ). (21) 

j>Q J 

This agrees with Eq. [lO] Here we do get the anticipated result. Note 
that this expression gives also a non-zero response for a perturbation 
which is non-zero only for a finite time, as opposed to Eq. |16| 

This sampling is the natural one fore deducing the general linear 
response theory. Doing the calculation for constant forcing does not 
elucidate the relevance of the choice of sampling. This sampling cor- 
responds to a Gcdankcnexperiment where the system is prepared in 
the distant past and we observe the difference of the perturbed and 
unperturbed evolution up to a given instant T. 
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3 Numerics 



3.1 Reick's formula 

For perturbations that are separable (X(t,x) = <fi(t)x(x)) and have a 
single driving frequency SI (<j>(t) — ecos(Qt)), the following sampling 
scheme for computing the susceptibility for a given observable A has 
been proposed by Reick [18] : 

1 N r 

k A (n) = l^lim— 5> <nt / p(dx) (AiftW-Atfix))) . (22) 
e °° t=i ^ 

This formula has been later adopted to analyze the output of a simple 
climate model |16j and a generalization has been proposed to study the 
nonlinear susceptibilities describing harmonic generation [14] . Appli- 
cability of this formula depends on performing numerical experiments 
where the initial samples approximate the unperturbed SRB measure 
P- 

Using our previous calculations, we want to circumstantiate the va- 
lidity of the formula. We apply Ruelle's response theory to obtain a 
perturbative expression of Eq. [22] in terms of quantities of the unper- 
turbed dynamics. 

1 N r 
l im Jj E ^ / ( A (fo(*)) ~ A(/*(»))J 

N t 

« JToo n E e ^ / E(/ J ')*« x 0'> f 3 (*)) D ( A ° ■ 

t=l J 3=1 



Here we encounter the same problem as in Eq. |16[ namely the coupling 
of the averages of the measure and the perturbation. Indeed, using 
Reick's formula sampling from an initial measure different from the 
unperturbed SRB measure, one does not get a reasonable response, as 
reported in [IB] , 

By sampling according to the unperturbed SRB measure p instead 
of I, the above equation becomes 



1 N N-j 

j™, Jj E E e< °* e<tW / P^W, X ) D ( A /*)(*) 

N 

= lim - / P {dx)G A {u, x) V e^X(j, z) . (23) 



3=1 

We insert the inverse discrete time Fourier transform 
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into Eq. [23] 

n 



/ P (dx)G A (oj, x) lim 1 ]T e in ^(i, i) 



JY 



i^oo N 4-* 2tt J_„ 

.7 = 1 



p{dx)G A {u,x) lim — > e lUj — / xSe^'du 

3- 



= I p{dx)G a{u : x) lim — f X(u>, x)un(£1 — ui)dui 
n->oo 2tt J_ 7T 

i r 

lira — / du!Uf<f(Q — oj)ka{u>) ■ (24) 



where 



u N (n -lS) = — ^ 
This can be rewritten by making use of 



N 

N ■ ' 



at la; 
•a; = x 



AT 



1-X 

as 

v AT 1 - e l (°-") 

which converges to as TV goes to infinity, except for Q = w. At f2 = o> 
the sum over j gives N. Hence, if X(ui,x) is integrable, we can take 
the limit in Eq. 24 inside the integral 

i r - 

— J X(uj,x)l {n} (uj)duj = 0, 

where is the indicator function on {fi}. We deduce that in the case 
of a general perturbation with a continuous Fourier spectrum, Reick's 
numerical approach cannot be applied. Note also that for finite time 
steps N (as is always the case for numerical experiments), there is an 
additional broadening of the signal of order 1 /N, as is apparent from 
Eq.[24] 

If on the other hand the Fourier transform is singular, for example 
X(uj,x)=S{n-cj) X (x), 
which corresponds to the monochromatic signal 

we have that 

N 



lim lye^XU,x) = ^ 
jv^oo x ^ u ' ; 2tt 

3 = 1 
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Therefore, Eq. [23] becomes 

J p{dx)G A (u,x)x{x) = «a(^) j 

as predicted by Reick. 

The above calculation shows how the explicit expansion of Reick's 
response formula allows us to interpret its finite time behaviour. Eq. |24| 
shows how this sampling scheme amounts to filtering the susceptibility 
ka with the function u^. 

Note that in the case of a several frequencies contributing to the 
forcing, we are again in the case of general periodic forcing. The sus- 
ceptibility can in this case be computed in two ways. Either one uses 
Reick's formula at every frequency present in the signal, which amounts 
to doing spectroscopy. On the other hand, one can also compute the 
full response StP at all phases over on period and apply a Fourier trans- 
form to this time-dependent function. The response at any one specific 
phase can be efficiently computed with the periodic sampling strategy 
proposed in Eq. [18] In this approach each value for the difference of A 
between perturbed and unperturbed is processed only once, compared 
to the summation being done for every frequency with Reick's formula. 

3.2 Sampling continuous spectra 

The discussion in the previous subsection demonstrates how sampling 
according to Formula [22] can only give a correct result in cases where 
the Fourier spectrum of the perturbation is discrete. In this subsection 
we explore how the discussion on the expansion at a fixed time hori- 
zon can help us find a sampling for the case of a continuous Fourier 
spectrum. 

One possibility is to sample the response directly from the full for- 
mula of the perturbation of the SRB measure: 

6rp(A) = lim \ V / l(dx)A(f T 0...0 f T . k {x)) - A(f k (x)) 

k=l J 

This sampling however entails some practical difficulties. As can be 
seen from the formula, an ergodic average is taken over the length of 
the numerical run k. Increasing k to k + 1 is equivalent to altering 
the initial conditions from x to /r-fe-i(x). This trajectory cannot 
be recovered from the previously calculated trajectories of length k. 
Hence one needs to redo the calculations of the trajectories for every 
value of k. As we will see, less costly sampling methods can be devised. 

To study the behaviour of different sampling methods, let us define 
the following quantity: 

5^p T {A) = J l(dx)A(f T o . . . o f T _ k o p{x)) - A(f k+n (x)) 

By changing k, we control the length of time over which we observe 
the difference between the perturbed and unperturbed dynamics. The 
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initial measure I is furthermore transformed by n applications of the 
unperturbed dynamics. By increasing n, the initial measure I converges 
to the unperturbed SRB m easu re p. 



Making use of Equation 20 we can expand S^' 71 ^ p T to get a better 
idea of the behaviour of this quantity under different limits and ergodic 
averages: 

k-l r 

8^p T {A)=Y J / f {k ~ j+n) *l{dx)X{T - j,x)D{Ao f){x) (25) 

3=0 J 

The aim when constructing a sampling scheme is to take limits and 
ergodic averages over k and n in such a way that the response given 



by Equation 10 is obtained. As we have seen with Reick's formula, 
this convergence can depend on the perturbation X. Furthermore, as 
exemplified by the discussion in this section, numerical cost should be 
considered. The following ergodic mean: 

Jim - Y^5^pr{A) (26) 

t-*oo t *■ — ' 
k = l 



converges to the response 5tp(A) given by Eq. 25 for any value of n. 
This can be shown following the discussion provided in Section |2.2.2[ 
where we discuss the case n = 0. Another possibility to obtain the 



SRB measure p in Eq. 25 is to take the limit of n going to infinity for 
a fixed k. We obtain: 

fe-i . 

lim 5 {k ' n) PT (A) = V / p(dx)X(T - j, x)D(A o f)(x). (27) 

3=0 J 

where we have assumed that linim^oo f m *l = p. This expression tends 
to Stp(A) in the limit of k oo. From a theoretical point of view, an 
increase in the value of n simply translates into a change in the initial 
measure I. Numerically, though, doing a long initial unperturbed run 
will evolve the initial measure towards the invariant measure p, hence 



improving convergence when Eqs. (25l-(27) are considered. In fact, 
there are a number of different options when attempting to reach a 
good numerical convergence. These include 

• increasing the length of unperturbed and perturbed trajectories 
(n and k) 

• enlarging the number of initial conditions (chosen according to /) 

• deciding whether or not ergodic averaging is performed over n 
and k. 

In the limit of infinitely long perturbed runs, these approaches give the 
same result. However, for finite time they will perform differently. 



4 Conclusions 

In this paper we have reconsidered Ruelle's linear response theory by 
analyzing the impact of choosing different methods of sampling in re- 
lation to different classes of forcings. Explicitly doing an expansion of 
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the perturbed dynamics around the unperturbed measures allows us to 
explore which sampling methods converge and under which conditions. 

The general response formula is obtained by choosing a specific 
sampling where the system is prepared in the distant past and we ob- 
serve the difference of the perturbed and the unperturbed dynamics up 
to a given time T. By proposing a general decomposition of space-time 
dependent forcings using a Schauder decomposition, we have elucidated 
that it is possible to define elementary linear susceptibilities that allow 
to construct the response of the system to any pattern of forcing. 

The other possible sampling strategy, where the time horizon is 
not fixed, does not give rise to a natural response theory except for 
constant perturbations. In the case of periodic forcings one can obtain 
a meaningful formula by redefining appropriately the response, finely 
tuned to the forcing under investigation. One needs to subsample the 
signal with the same period of the forcing and explore all the initial 
phases. By taking this approach, it is in principle possible to discover 
the fundamental period of the external perturbations by varying the 
sampling period. Thanks to our approach we get a deeper understand- 
ing of the range of applicability of Reick's formula, which has been used 
as a signal processing tool to study the linear response of numerical 
models. 

Nonetheless, this approach fails if the forcing is not periodic, in 
which case we must resort to the fixed-time horizon framework to get 
a meaningful answer. In fact, our findings explain why considering 
the fixed-time horizon it is possible to analyze a response to forcings 
that have a continuous Fourier spectrum. The clarifications presented 
in this paper may be of relevance for devising the data processing for 
actual laboratory experiments on nonlinear systems. 

We also clarify that it is crucial in practical terms to use an ensem- 
ble approach where the initial conditions sample approximately the 
unperturbed SRB measure. Our calculation is explicitly performed 
for discrete time, but the analogous results for continuous time are 
presented in Appendix [XJ Moreover our considerations seem to be ap- 
propriate also for the case of nonlinear response IS] • 

To summarize, we have shown the following 

• Sampling a general response from an initial time up to a moving 
time horizon does not lead to a well-defined sampling method. 

• Starting the simulation at times in the distant past and averaging 
the response at a fixed time horizon always results in the full 
response of the system at the fixed point in time. This approach 
can be computationally inefficient. 

• Sampling a periodic response with a moving time horizon results 
in a response of the system as if it were forced with an averaged 
forcing. 

• In the periodic case, the full response can be computed by sam- 
pling with a horizon moving forward in time with steps of one 
period. This response depends on the initial phase. The sus- 
ceptibility can be computed through a Fourier transform of the 
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response. 

• A constant forcing can be considered as a periodic forcing with 
period 1 and can thus be sampled with a horizon moving with 
time steps of 1. 

• For periodic forcings, Reick's spectroscopic formula also allows 
to discern the response at different frequencies, i.e. the suscepti- 
bility. It gives a zero susceptibility for forcings with a continuous 
spectrum. 

We believe that the results presented in this article can be of interest 
to researchers interested in studying the response of complex systems 
to modulations of their internal parameters or to external perturba- 
tions. For various reasons, climate science is an especially promising 
field of application. First of all, we clarify crucial differences between 
sampling periodic and aperiodic forcings. This is a crucial issue if one 
wants to apply linear response theory to study different scenarios such 
as the response of the system to monotonically increasing CO2 levels 
(see a forthcoming paper by the authors) versus its response to periodic 
forcings such as those due to astronomical and astrophysical phenom- 
ena. In particular, we have proposed a parsimonious but effective way 
for analysing periodic - but non-monochromatic - forcings. Again, this 
setting is applicable to climate science due to the presence of cycles 
with different time scale, such as the daily, yearly and solar cycles. 

Future work will address the investigation of the response of a non- 
equilibrium system to a general random field. Moreover, we will ana- 
lyze the impact of the various sampling schemes described in this paper 
when studying the output of numerical models. 
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A Continuous time response formulas 

Here we give the formulas for continuous time systems corresponding 
to the ones presented in the main text. The time evolution is in this 
setting given by a differential equation 

dx cv \ 
* = F{X) ' 

resulting in a flow x(t + s) = f s (x(t)). The SRB measure is given by 

p(A) = lim 1 / da f l(dx)A(f s {x)) . 
t Jo J 

For a separable perturbation 

dx 

— =F{x)+ X {x)4>{t) 
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the susceptibility is given by 

poo p 

- J dte 1 ^ J p(dx) X (x)D(A o f)(x) . 

In case of a general perturbation F(x) — > F(x) + ir) the linear 
response becomes: 



<5 T p(A) ss jf dr J p{dx)X(T - r, x)D(A o / 1 
The SRB measure with a moving time horizon is: 



1 '* 



)(*) 



p T = lim - / ds / /(dx)A(/Z' +t (a;)) 
t-s-oo f ,/ J 

and the SRB measure with a fixed time horizon: 

p T = lim 7 / da / l(dx)A(ff_ t (x)) 

where /' 2 (a;) is a trajectory of the perturbed system, starting at time 
ti in x and evolving up to time t^. 
Reick's formula now becomes: 

k a (oj) = lim lim - f dte mt f p(dx) U(%(x)) - A{f*(x)j) 
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